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Abstract 

We present a family of numerical implementations of Kato's ODE propagating global 
bases of analytically varying invariant subspaces, of which the first-order version is a 
surprising simple "greedy algorithm" that is both stable and easy to program and 
the second-order version a relaxation of a first-order scheme of Brin and Zumbrun. 
The method has application to numerical Evans function computations used to assess 
stability of traveling-wave solutions of time-evolutionary PDE. 

1 Introduction 

Let P(A) G C nxn be a projection, P 2 = P, and subspace S(X) C C n its range, with P 
depending analytically on A within a simply connected subset A of the complex plane. 
Then, a standard result in matrix perturbation theory [K] is that there exists a global 
analytic basis {^(A)} of S on A; moreover, expressing r* as column vectors, this can be 
prescribed constructively as the solution R = (r±, . . . , r^) of Kato's ODE 

(1.1) R' = (P'P - PP')R, R(X ) = R°, 

where the initializing value R° is any matrix whose columns form a basis for S(Xo), and '"" 
denotes differentiation with respect to A. This prescription is also "minimal" in the sense 
that PR' = 0, i.e., the derivative of basis R lies entirely in the direction complementary to 
its span (the kernel of P); see [HSZ} IHLZ] for further discussion. 

The problem of computing such an analytically varying basis is important in numerical 
Evans function computations for determining stability of traveling waves. Roughly speaking, 
analytic bases for stable and unstable manifolds of certain limiting coefficient matrices are 
used to define an analytic Evans function, whose winding number around a contour T C A 
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counts the number of unstable eigenvalues enclosed by T of the linearized operator about 
the wave, with zero winding number corresponding to stability. For further discussion, 
see [GZl 113 |B?Zl IBDU1 iHgZl IHuZl IBHKZl iHLZl ICHNZl IHLyZlj |HLyZ2| and references 
therein. 

Various algorithms for numerical determination of bases R have been introduced in 
[BrZl IBDG} IHSZj . each of which turn out to be equivalent to (a discretization of) (jl.ip . and 
each of which is 0(n 3 ) in complexity (though with different coefficients). The purpose of 
this brief note is to introduce a new and particularly simple discretization of (jl.ip . which 
at the first order of accuracy consists of what might be called a "local greedy algorithm". 
Namely, choosing a set of mesh points Xj around T and denoting by Rj the approximation 
of R(Xj), our first-order scheme is simply 

(1.2) Rj+i = Pj+lRji Ro = R°- 

That is, the continuation of basis R to each new step is obtained by projecting the value at 
the previous step onto the new subspace: the simplest possible choice, and one that at first 
sight seems entirely local. However, remarkably, this local choice leads to a globally defined 
basis; in particular, upon traversing the entire contour T and returning to Xl = Ao, we find 
that, up to convergence error, the value of Rl returns to the starting value Rq = R°. This 
is both simpler and faster than its closest relative in jBrZ]; indeed, it is completely trivial 
to program (the main issue in most Evans function applications). 
This simplification is based on the reduced version 

(1.3) R' = P'R, R(X ) = R°, 

of (jl.ip ([K], pp. 99-101), which readily yields minimal difference schemes to all orders 
of accuracy. A particularly attractive version when speed is an issue appears to be the 
second-order version, which is a relaxation of the first-order scheme in |BrZ| . 



1.1 The reduced ODE 

We begin by recalling the properties of the reduced Kato ODE (|1.3j) . 

Proposition 1.1. There exists a global solution R of (jl.3p on any simply connected domain 
A containing Ao, withrankR(X) = rankR . Moverover, if P(Xq)R° = R°, then: (i) PR = R. 
(ii) PR' = 0. (in) R satisfies (fLTjl . 

Proof. As a linear ODE with analytic coefficients, (jl.3p possesses an analytic solution in 
a neighborhood of Ao, that may be extended globally along any curve, whence, by the 
principle of analytic continuation, it possesses a global analytic solution on any simply 
connected domain containing Ao [K]. Constancy of ranki?(A) follows likewise by the fact 
that R satisfies a linear ODE. 

Differentiating the identity P 2 = P following [K] yields PP' + P'P = P' , whence, 
multiplying on the right by P, we find the key property 

(1.4) PP'P = 0. 
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From (|1.4p . we obtain 

(PR - R)' = (P'R + PR' - R 1 ) = P'R + (P - I) P'R = PP'R, 

which, by PP'P = and P 2 = P gives 

(PR - R)' = -PP'(PR - R), (PR - R)(X ) = 0, 

from which (i) follows by uniqueness of solutions of linear ODE. Expanding PR' = PP'R 
and using PR = R and PP'P = 0, we obtain PR' = PP'PR = 0, verifying (ii). Finally, 
using (i) and (ii), we obtain R' = P'R = P'PR - PR' = (P'P - PP')R, verifying (iii). □ 

Remark 1.2. Conversely, (i)— (ii) imply (JOD through R' = (PR)' = P'R + PR' = P'R. 

1.2 Numerical implementation 
1.2.1 First-order version 

Approximating P'(Xj) to first order by the finite difference (Pj+i — Pj)/(Xj+i — Xj) and 
substituting this into a first-order Euler scheme gives 

Rj+i = Rj + (Xj+i - Xj)- 2 ^ -tRj, 

or Rj+i = Rj + Pj+iRj — PjRj, yielding greedy algorithm (jl.2p by the property PjRj = Rj 
(Note: preserved exactly by the scheme). 

Remark 1.3. The same procedure applied to the original equation (jl.ip yields 

R j+1 = Rj + (Pj+iPj - PjPj+i)Rj, 

or, following with a projection Pj+i to stabilize the scheme without changing the order of 
accuracy, the first-order scheme 

(1.5) Rj+i = Pj+lRj + Pj+l(Pj+lPj ~ PjPj+i)Pj = Pj+i[l + Pj(.I ~ -fj+i)) -Re- 
introduced in |BrZ| . This is slightly more costly at two evaluations of P on the average 
and three matrix multiplications vs. one evaluation of P and one matrix multiplication for 
(jl.2p . Depending on the cost of evaluating P, and whether or not the mesh is fixed (in 
which case evaluations of P may be shared), this can vary between approximately one and 
three times the cost of (jl.2p . 



3 



1.2.2 Second-order version 



To obtain a second-order discretization of (|1.3p . we approximate Rj + \—Rj « AXjP^ +1 ^ 2 Rj +1 / 2 , 
good to second order, where AXj := — Xj. Noting that Rj+1/2 ~ Pj+ifePj to second 
order, by (jl,2p . and approximating Pj+1/2 ~ |(Pj+i + Pj), also good to second order, and 
Pj + U2 ~ (-Pj'+l ~~ Pj)/AXj, we obtain, putting everything together and rearranging, 
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Rj+i = R, + 2(^+1 " p i)( p i+i + p j)«. 
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Stabilizing by following with a projection Pj+i, we obtain after some rearrangement the 
reduced second-order explicit scheme 

(1.6) R j+1 = P j+1 [I + ±P j (I-Pj +1 )]Rj, 

which may be recognized as a relaxation of the first-order scheme (|1.5p . 

This has the same computational cost as (jl.5|) . i.e., two evaluations of Pj and three 
matrix multiplications, while the number of steps goes as erance vs. 1 /Tolerance 

for first-order, so 10 times fewer for typical tolerance 10 -2 , for computational savings of ten 
times over (|1.5p and four times over (jl.2p . in the worst case (P inexpensive compared to 
matrix multiplication) that (II. 6p is three times as expensive as (11.2p . This is the version 
we recommend for serious computations. For individual numerical experiments the simpler 
greedy algorithm (11. 2p will often suffice (see discussion, Section fT~4l) . 

1.3 Third and higher-order versions 

Arbitrarily higher-order schemes may be obtained by Richardson extrapolation starting 
from scheme (ll.2p or (jl.6p . For example, second-order Richardson extrapolation applied to 
(|1.2p yields an alternative second-order scheme 

(1.7) R 3+l = P J+1 {2P ]+l/2 -I)R v 

while third-order extrapolation applied to (jl.7|) yields the third-order scheme 



(1.8) R j+1 = P 3 
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Third-order extrapolation applied to (|1.6p yields the simpler third-order scheme 
(1.9) 



Rj. 



Here, fractional indices denote points along the line segment between Xj and Xj + i/ 2 with 
corresponding fractional distance. 
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More generally, denoting by T m,h the matrix advancing Rj to Rj+i for an mth order 
scheme with step h £ C, we obtain by Richardson extrapolation an (m + l)st order scheme 

j ~ 2 m — 1 J +1 J ~~ 2 m — 1 J 

When P is costly, it appears preferable to use schemes involving evaluations at only integer 
steps, in order to share evaluations of P (of course, this assumes a fixed, or non-adaptive, 
mesh, which may or may not be desirable). For example, explicit mth-order Euler approxi- 
mating derivatives of P at Xj by Lagrange interpolation yields an (m + l)-step scheme with 
integer steps. However, the complexity of the resulting schemes makes these unappealing 
in practice. 

Among higher-order schemes, we thus suggest only (|1.6p . or for strict tolerance (jl.9p . 



1.4 Implementation in numerical Evans function computations 

For Evans computations, P(X) is the eigenprojection onto the stable (unstable) subspace of 
a given matrix A(X). Thus, it may be prescribed uniquely as 

(1.10) P(X) = R(L* R)^ 1 L* , 

for any choice of right and left bases R and L (matrices whose columns consist of basis 
elements, as before). 

Ordered Schur decomposition, an 0(n 3 ) operation supported as an automatic function 
in programming packages such as MATLAB, applied to A and A*, respectively, gives or- 
thonormal right and left bases of the left and right stable subspaces of A, hence an optimally 
conditioned choice in (ll.lOp . Thus, evaluation of P is in practice straightforward to pro- 
gram. On the other hand, it is typically an expensive (i.e., large coefficient) 0(n 3 ) call 
involving Schur decomposition and several matrix multiplications, so that it is desirable 
to minimize the number of evaluations of P in numerical continuation algorithms. For a 
fixed, i.e., non-adaptive, mesh (at least for (II. 2ft . (II. 6j) . or explicit Euler schemes that are 
evaluated at mesh points only), P need be evaluated only once for every mesh point, so 
that higher-order schemes are clearly preferable in this application. On the other hand, 
evaluation of the Evans function, involving solution of a further variable-coefficient ODE 
initialized with R, costs far more than the computation of R, so that these details can be 
ignored in most computations with (relatively) little effect. Ordered Schur decomposition 
(MATLAB version) has been used with good results in [HuZl HLyZ2| . 
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